! INLET/OUTLET boundary conditions:
! Total temperature and total pressure are given
! VSemiletov, 2012-01-10

! inlet
subroutine CalcInletBoundCond_TotalPresTempGiven( Un,Ucell,Pn,Rhon,meshSf,Pin,Tin,Uin)
use DataPerfectGas, only: cv, Gamma
implicit none

    real(8),intent(in) :: Ucell(3),meshSf(3),Pin,Tin,Uin
    real(8),intent(out) :: Un(3),Pn,Rhon

    real(8) :: c,M,Ptotal,Ttotal,Tn, Un2

    c = sqrt(Gamma*(Gamma-1d0)*cv*Tin) ! Pin/RhoIn)
    M = Uin / c

    Ptotal = Pin * (1d0 + 0.5d0*(Gamma-1d0)*M**2)**(Gamma/(Gamma-1d0))
    Ttotal = Tin * (1d0 + 0.5d0*(Gamma-1d0)*M**2)

    Un(1:3) = Ucell(1:3)
    Un2 = sum(Un(1:3)*Un(1:3)) !meshSf(1:3))
!    Un(1:3) = sqrt(Un2)*meshSf(1:3)

    Tn = Ttotal - 0.5d0*Un2/(Gamma*cv)
    Pn = Ptotal / (1d0 + 0.5d0*Un2/(Gamma*(Gamma-1d0)*cv*Tn))**(Gamma/(Gamma-1d0))
    Rhon = Pn / ((Gamma-1d0)*cv*Tn)

!    un(1) = ucsb(1)
!    Tn = Ttotal - 0.5d0*un(1)**2/(Gamma*cv)
!    Pn = Ptotal / (1d0 + 0.5d0*un(1)**2/(Gamma*cv*Tn))**(Gamma/(Gamma-1d0))
!    Ron = Pn / ((Gamma-1d0)*cv*Tn)

    return
end subroutine CalcInletBoundCond_TotalPresTempGiven

! outlet
subroutine CalcOutletBoundCond_TotalPresTempGiven(Un,Presn,Rhon,Ucell,PresCell,RhoCell,Uout)
use DataPerfectGas !, only: cv, Gamma
implicit none

    real(8),intent(in) :: PresCell,RhoCell,Uout,Ucell(3)
    real(8),intent(out) :: Un(3), Presn, Rhon
    
    real(8) :: Tcell, c2,M2,Ucell2, Tn, Ptotal,Ttotal,Un2

    Un(1) = Uout
    Un(2:3) = 0d0
    Presn = PresCell
    Rhon = RhoCell
    
    Tcell = TemperaturePRho(PresCell,RhoCell)

    c2 = Gamma*(Gamma-1d0)*cv*Tcell ! Pin/RhoIn)
    Ucell2 = sum(Ucell**2)
    M2 = Ucell2 / c2

    Ptotal = PresCell * (1d0 + 0.5d0*(Gamma-1d0)*M2)**(Gamma/(Gamma-1d0))
    Ttotal = Tcell * (1d0 + 0.5d0*(Gamma-1d0)*M2)

!    Un(1:3) = Ucell(1:3)
    Un2 = sum(Un(1:3)*Un(1:3)) !meshSf(1:3))
!    Un(1:3) = sqrt(Un2)*meshSf(1:3)

    Tn = Ttotal - 0.5d0*Un2/(Gamma*cv)
    Presn = Ptotal / (1d0 + 0.5d0*Un2/(Gamma*(Gamma-1d0)*cv*Tn))**(Gamma/(Gamma-1d0))
    Rhon = Presn / ((Gamma-1d0)*cv*Tn)

    return
end subroutine CalcOutletBoundCond_TotalPresTempGiven